pacman::p_load("adehabitatHR","data.table","ggfortify","OpenStreetMap","grid","move","moveVis","OpenStreetMap","pbapply","maptools","rgdal","plotly","sp","tidyverse","viridis")

1 Import Data Set

This is the same data set used for the Basic Mapping assignment. These GPS points were collected for three Louisiana Waterthrush territories at Cheatham Wildlife Management Area in Cheatham Co, TN. For the purpose of getting a later line of code to work, I changed the ID’s for each bird to LOWA1-LOWA3.

data <- read.csv("./Data/lowa_location_data.csv")
head(data)

2 QAQC Plot

qaqc_plot <- ggplot() + geom_point(data=data, 
                                   aes(Easting,Northing,
                                       color=ï..id)) +
                        labs(x="Easting", y="Northing") +
                        guides(color=guide_legend("Identifier"))

ggplotly(qaqc_plot)
lapply(split(data, data$ï..id), 
       function(x)write.csv(x, file = paste(x$ï..id[1],".csv", sep = ""), row.names = FALSE))
files <- list.files(path = ".", pattern = "[LOWA]+[0-9]+", full.names = TRUE)
utm_points <- cbind(data$Easting, data$Northing)
utm_locations <- SpatialPoints(utm_points, 
                 proj4string=CRS("+proj=utm +zone=16 +datum=WGS84"))
proj_lat.lon <- as.data.frame(spTransform(
                utm_locations, CRS("+proj=longlat +datum=WGS84")))
colnames(proj_lat.lon) <- c("x","y")
raster <- openmap(c(max(proj_lat.lon$y)+0.01, min(proj_lat.lon$x)-0.01), 
                  c(min(proj_lat.lon$y)-0.01, max(proj_lat.lon$x)+0.01), 
                  type = "bing")
raster_utm <- openproj(raster, 
              projection = "+proj=utm +zone=16 +datum=WGS84 +units=m +no_defs")

\(~\)

3 Plot the Data

autoplot.OpenStreetMap(raster_utm, expand = TRUE) + theme_bw() +
  theme(legend.position="bottom") +
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_point(data=data, aes(Easting,Northing,
             color=ï..id), size = 3, alpha = 0.8) +
  theme(axis.title = element_text(face="bold")) + labs(x="Easting",
        y="Northing") + guides(color=guide_legend("Identifier"))

\(~\)

4 Minimum Convex Plygon

mcp_raster <- function(filename){
  data <- read.csv(file = filename)
  x <- as.data.frame(data$Easting)
  y <- as.data.frame(data$Northing)
  xy <- c(x,y)
  data.proj <- SpatialPointsDataFrame(xy,data, proj4string = CRS("+proj=utm +zone=16 +datum=WGS84 +units=m +no_defs"))
  xy <- SpatialPoints(data.proj@coords)
  mcp.out <- mcp(xy, percent=100, unout="ha")
  mcp.points <- cbind((data.frame(xy)),data$ï..id)
  colnames(mcp.points) <- c("x","y", "identifier")
  mcp.poly <- fortify(mcp.out, region = "id")
  units <- grid.text(paste(round(mcp.out@data$area,2),"ha"), x=0.85,  y=0.95,
                     gp=gpar(fontface=4, col="white", cex=0.9), draw = FALSE)
  mcp.plot <- autoplot.OpenStreetMap(raster_utm, expand = TRUE) + theme_bw() + theme(legend.position="none") +
    theme(panel.border = element_rect(colour = "black", fill=NA, size=1)) +
    geom_polygon(data=mcp.poly, aes(x=mcp.poly$long, y=mcp.poly$lat), alpha = 0.5, fill = "red") +
    geom_point(data=mcp.points, aes(x=x, y=y)) + 
    labs(x="Easting (m)", y="Northing (m)", title=mcp.points$identifier) +
    theme(legend.position="none", plot.title = element_text(face = "bold", hjust = 0.5)) + 
    annotation_custom(units)
  mcp.plot
}

pblapply(files, mcp_raster)
## [[1]]

## 
## [[2]]

## 
## [[3]]

\(~\)

5 Kernel Density Estimate

kde_raster <- function(filename){
  data <- read.csv(file = filename)
  x <- as.data.frame(data$Easting)
  y <- as.data.frame(data$Northing)
  xy <- c(x,y)
  data.proj <- SpatialPointsDataFrame(xy,data, proj4string = CRS("+proj=utm +zone=16 +datum=WGS84 +units=m +no_defs"))
  xy <- SpatialPoints(data.proj@coords)
  kde<-kernelUD(xy, h="href", kern="bivnorm", grid=100)
  ver <- getverticeshr(kde, 95)
  kde.points <- cbind((data.frame(data.proj@coords)),data$ï..id)
  colnames(kde.points) <- c("x","y","identifier")
  kde.poly <- fortify(ver, region = "id")
  units <- grid.text(paste(round(ver$area,2)," ha"), x=0.85,  y=0.95,
                     gp=gpar(fontface=4, col="white", cex=0.9), draw = FALSE)
  kde.plot <- autoplot.OpenStreetMap(raster_utm, expand = TRUE) + theme_bw() + theme(legend.position="none") +
    theme(panel.border = element_rect(colour = "black", fill=NA, size=1)) +
    geom_polygon(data=kde.poly, aes(x=kde.poly$long, y=kde.poly$lat), alpha = 0.5, fill = "red") +
    geom_point(data=kde.points, aes(x=x, y=y)) +
    labs(x="Easting (m)", y="Northing (m)", title=kde.points$identifier) +
    theme(legend.position="none", plot.title = element_text(face = "bold", hjust = 0.5)) + 
    annotation_custom(units)
  kde.plot
}

pblapply(files, kde_raster)
## [[1]]

## 
## [[2]]

## 
## [[3]]

The three KDE maps have large buffers and aren’t very accurate, mostly due to the low number of data points to work from.